setwd("E://OneDrive//研究//3-小实验TU//5-论文//入侵程度//数据与分析")

library(xlsx)
library(tidyr)
library(plyr)
library(dplyr)
library(car)
library(agricolae)
library(ggplot2)
library(ggpubr)
library(gridExtra)

##################
###  数据输入  ###
##################
rm(list=ls())

grow <- read.xlsx("20221024.xlsx", sheetName = "grow" , encoding = "UTF-8")
grow$Treat <- as.factor(grow$Treat)
grow$T <- as.factor(grow$T)
grow$U <- as.factor(grow$U)
grow$ID <- as.factor(grow$ID)
grow$PT <- as.factor(grow$PT)

leaf <- read.xlsx("20221024.xlsx", sheetName = "leaf" , encoding = "UTF-8")
leaf$Treat <- as.factor(leaf$Treat)
leaf$T <- as.factor(leaf$T)
leaf$U <- as.factor(leaf$U)
leaf$ID <- as.factor(leaf$ID)
leaf$PT <- as.factor(leaf$PT)

ps <- read.xlsx("20221024.xlsx", sheetName = "photosynthesis" , encoding = "UTF-8")
ps$Treat <- as.factor(ps$Treat)
ps$T <- as.factor(ps$T)
ps$U <- as.factor(ps$U)
ps$ID <- as.factor(ps$ID)
ps$PT <- as.factor(ps$PT)

#soil <- read.xlsx("20221024.xlsx", sheetName = "soil" , encoding = "UTF-8")
#soil$Treat <- as.factor(soil$Treat)
#soil$T <- as.factor(soil$T)
#soil$U <- as.factor(soil$U)
#soil$ID <- as.factor(soil$ID)

##################
###  定义函数  ###
##################
Normal_t.yb <- function(data, Col_DV, Col_IV, Col_group = NULL) {
  library(psych)
  library(plyr)
  library(car)
  df <- data
  if (!is.null(Col_group)) {
    df <- unite(df, "group", all_of(Col_group), sep = "_", remove = F)
    df$group <- as.factor(df$group)
    group <- unique(df$group)
  } else {
    df$group <- as.factor(1)
    group <- unique(df$group)
  }
  df <- unite(df, "IV", all_of(Col_IV), sep = "_", remove = F)
  df$IV <- paste("IV", df$IV, sep = "_")
  df$IV <- as.factor(df$IV)
  out <- list()
  out[[1]] <- df
  out[[2]] <- data.frame(group = "", DV = "", Conversion_process = "")[-1,]
  for (i in c(1:length(group))) {
    df1 <- df[df$group == group[i],]
    df1$IV <- as.factor(df1$IV)
    IV <- unique(df1$IV)
    for (j in c(1:length(Col_DV))) {
      homogeneity <- c()
      for (r in 1:4) {
        if (r == 1)
          df1$y <- df1[,Col_DV[j]]
        if (r == 2)
          df1$y <- df1[,Col_DV[j]]^0.5
        if (r == 3)
          df1$y <- log(df1[,Col_DV[j]])
        if (r == 4)
          df1$y <- 1/df1[,Col_DV[j]]
        homogeneity[r] <- as.data.frame(leveneTest(eval(parse(text = paste("y", paste(Col_IV, collapse = "*"), sep = "~"))), data = df1, center = mean))[1,3]
      } 
      if (homogeneity[1] > 0.05)
        Homogeneity_variance <- "Homogeneity"
      else {
        if (any(homogeneity[2:4] > 0.05)) {
          r <- which(homogeneity[2:4] == max(homogeneity[2:4]))+1
          if (r == 2)  {
            df1[,Col_DV[j]] <- df1[,Col_DV[j]]^0.5
            Homogeneity_variance <- "+SQRT"
          }
          if (r == 3) {
            df1[,Col_DV[j]] <- log(df1[,Col_DV[j]])
            Homogeneity_variance <- "+Ln"
          }
          if (r == 4) {
            df1[,Col_DV[j]] <- 1/df1[,Col_DV[j]]
            Homogeneity_variance <- "+1/x"
          }
        } else {
          Homogeneity_variance <- "Non homogeneity"
        }
      }
      Is_normal <- F
      r <- 0
      while (all((!Is_normal), r < 10)) {
        r <- r + 1
        normal <- data.frame(group = "", DV = "", IV = "", n = "", KS = "", KS_P = "", SW = "", SW_P = "")[-1,]
        for (k in 1:length(IV)) {
          sample <- df1[df1$IV == IV[k], Col_DV[j]]
          n <- length(sample)
          ks <- ks.test(sample, pnorm, mean(sample), sd(sample))
          sw <- shapiro.test(sample)
          normal <- rbind(normal, c(as.character(group[i]), as.character(Col_DV[j]), as.character(IV[k]), n, 
                                    round(ks[[1]], 4), round(ks[[2]], 4), 
                                    round(sw[[1]], 4), round(sw[[2]], 4)))
        }
        names(normal) <- c("group", "DV", "IV", "n", "KS", "KS_P", "SW", "SW_P")
        if (mean(as.numeric(normal$n)) < 100) {
          Is_normal <- all(normal$SW_P > 0.05)
        } else {
          Is_normal <- all(normal$KS_P > 0.05)
        }
        if (!Is_normal) {
          Skewness <- data.table::rbindlist(
            describe(eval(parse(text = paste(Col_DV[j], "IV", sep = "~"))), data = df1), idcol="IV")
          Skewness$skew.se <- sqrt(6 * Skewness$n * (Skewness$n - 1)/(Skewness$n - 2)/(Skewness$n + 1)/(Skewness$n + 3))
          Skewness$skew.Z <- Skewness$skew/Skewness$skew.se
          if (mean(Skewness$n) <= 50) 
            z_0.05 <- 1.96
          if (mean(Skewness$n) <= 300 & mean(Skewness$n) > 50) 
            z_0.05 <- 3.29
          if (mean(Skewness$n) > 300) 
            z_0.05 <- mean(2/Skewness$skew.se)
          if (abs(mean(Skewness$skew.Z)) < z_0.05) {
            if (r == 1)
              Conversion_process <- "skew ∈ normal"
            else
              Conversion_process <- paste(Conversion_process, " skew∈normal", sep = " ~")
            Is_normal <- T
          } else {
            Conversion_process <-""
            if (mean(Skewness$skew) > 0) {
              if (abs(mean(Skewness$skew.Z)) < 3*z_0.05)  {
                df1[,Col_DV[j]] <- df1[,Col_DV[j]]^0.5
                Conversion_process <- paste(Conversion_process, " +SQRT", sep = " ~")
              }
              if (abs(mean(Skewness$skew.Z)) >= 3*z_0.05 & abs(mean(Skewness$skew.Z)) < 7*z_0.05) {
                df1[,Col_DV[j]] <- log(df1[,Col_DV[j]])
                Conversion_process <- paste(Conversion_process, " +Ln", sep = " ~")
              }
              if (abs(mean(Skewness$skew.Z)) >= 7*z_0.05) {
                df1[,Col_DV[j]] <- 1/df1[,Col_DV[j]]
                Conversion_process <- paste(Conversion_process, " +1/x", sep = " ~")
              }
            }
            if (mean(Skewness$skew) < 0) {
              if (abs(mean(Skewness$skew.Z)) < 3*z_0.05)  {
                df1[,Col_DV[j]] <- (max(df1[,Col_DV[j]])+1-df1[,Col_DV[j]])^0.5
                Conversion_process <- paste(Conversion_process, " -SQRT", sep = " ~")
              }
              if (abs(mean(Skewness$skew.Z)) >= 3*z_0.05 & abs(mean(Skewness$skew.Z)) < 7*z_0.05) {
                df1[,Col_DV[j]] <- log(max(df1[,Col_DV[j]])+1-df1[,Col_DV[j]])
                Conversion_process <- paste(Conversion_process, " -Ln", sep = " ~")
              }
              if (abs(mean(Skewness$skew.Z)) >= 7*z_0.05) {
                df1[,Col_DV[j]] <- 1/(max(df1[,Col_DV[j]])+1-df1[,Col_DV[j]])
                Conversion_process <- paste(Conversion_process, " -1/x", sep = " ~")
              }
            }
          }
        } else {
          if (r == 1)
            Conversion_process <- "Normal"
          else
            Conversion_process <- paste(Conversion_process, " Normal", sep = " ~")
        }
      }
      if (r == 10)
        Conversion_process <- paste(Conversion_process, " 10times", sep = " ~")
      out[[1]][(out[[1]]$IV %in% df1$IV & out[[1]]$group == df1$group), Col_DV[j]] <- df1[,Col_DV[j]]
      out[[2]] <- rbind(out[[2]], c(group[i], Col_DV[j], Homogeneity_variance, Conversion_process))
    }
  }
  names(out[[2]]) <- c("group", "DV", "Homogeneity_variance", "Conversion_process")
  out[[1]] <- out[[1]][,colnames(data)]
  return(out)
}  
#正态性检验；输入(列名)：data-数据表(数据框)、DV-因变量(向量)、VI-自变量/固定因子(向量)、group-多重比较分组(向量)
#输入：out[[1]]：转换后的数据与输入的data相同；out[[2]]：转换过程(Homogeneity_variance=方差齐性；Conversion_process=正态性转换)

ANOVA.yb <- function(df, Col_DV, Col_IV, Col_group = NULL) {
  library(tidyr)
  library(car)
  library(plyr)
  library(dplyr)
  library(agricolae)
  
  if (!is.null(Col_group)) {
    df <- unite(df, "group", all_of(Col_group), sep = "_", remove = F)
    df$group <- as.factor(df$group)
    group <- unique(df$group)
  } else {
    df$group <- as.factor(1)
    group <- unique(df$group)
  }
  
  result <- list()
  
  for (i in c(1:length(group))) {
    result[[i]] <- list()
    df1 <- df[df$group == group[i],]
    df1 <- unite(df1, "IV", all_of(Col_IV), sep = "_", remove = F)
    df1$IV <- paste("IV", df1$IV, sep = "_")
    df1$IV <- as.factor(df1$IV)
    IV <- unique(df1$IV)
    
    for (j in c(1:length(Col_DV))) {
      result[[i]][[j]] <- list()
      length2 <- function (x, na.rm=T) {
        if (na.rm) sum(!is.na(x))
        else       length(x)
      }
      datac <- ddply(df1, c(Col_IV, "IV"), .drop=T,
                     .fun = function(xx, col) {
                       c(N    = length2 (xx[[col]], na.rm=T),
                         mean = mean    (xx[[col]], na.rm=T),
                         sd   = sd      (xx[[col]], na.rm=T),
                         min  = min     (xx[[col]], na.rm=T),
                         max  = max     (xx[[col]], na.rm=T)
                       )
                     },
                     Col_DV[j]
      )
      datac$se <- datac$sd / sqrt(datac$N) 
      datac$down <- datac$mean - datac$se
      datac$up <- datac$mean + datac$se
      result[[i]][[j]][[3]] <- datac
      
      model <- paste(Col_DV[j], paste(Col_IV, collapse = "*"), sep = "~")
      normal <- data.frame(group = "", DV = "", IV = "", KS = "", KS_P = "", SW = "", SW_P = "")[-1,]
      for (k in c(1:length(IV))) {
        sample <- df1[df1$IV == IV[k], Col_DV[j]]
        ks <- ks.test(sample, pnorm, mean(sample), sd(sample))
        sw <- shapiro.test(sample)
        normal <- rbind(normal, c(as.character(group[i]), as.character(Col_DV[j]), as.character(IV[k]), 
                                  round(ks[[1]], 4), round(ks[[2]], 4), 
                                  round(sw[[1]], 4), round(sw[[2]], 4)))
      }
      names(normal) <- c("group", "DV", "IV", "KS", "KS_P", "SW", "SW_P")
      result[[i]][[j]][[1]] <- as.data.frame(leveneTest(eval(parse(text = model)), data = df1, center = mean))
      result[[i]][[j]][[2]] <- as.data.frame(summary(aov(eval(parse(text = model)), data = df1))[[1]])
      SNK <- SNK.test(aov(eval(parse(text = paste(Col_DV[j], "IV", sep = "~"))), data = df1), "IV")$group
      SNK <- data.frame(IV = rownames(SNK), SNK = SNK$groups)
      LSD <- LSD.test(aov(eval(parse(text = paste(Col_DV[j], "IV", sep = "~"))), data = df1), "IV", p.adj="none")$group
      LSD <- data.frame(IV = rownames(LSD), LSD = LSD$groups)
      duncan <- duncan.test(aov(eval(parse(text = paste(Col_DV[j], "IV", sep = "~"))), data = df1), "IV")$group
      duncan <- data.frame(IV = rownames(duncan), duncan = duncan$groups)
      HSD <- HSD.test(aov(eval(parse(text = paste(Col_DV[j], "IV", sep = "~"))), data = df1), "IV")$group
      HSD <- data.frame(IV = rownames(HSD), HSD = HSD$groups)
      mc <- merge(merge(SNK, LSD, by = "IV"), merge(duncan, HSD, by = "IV"), by = "IV")
      result[[i]][[j]][[3]] <- merge(result[[i]][[j]][[3]], merge(normal, mc, by = "IV"), by = "IV")
      result[[i]][[j]][[3]] <- result[[i]][[j]][[3]][, c("group", "DV", Col_IV, "N", "mean", "sd", "min", "max", "se", "down", "up", "KS", "KS_P", "SW", "SW_P", "SNK", "LSD", "duncan" ,"HSD")]
    }
  }
  return(result)
} 
#方差分析；输入(列名)：df-数据表(数据框)、DV-因变量(向量)、VI-自变量/固定因子(向量)、group-多重比较分组(向量)
#输入：result[[group]][[DV]][[1-levene(n>50-KS); 2-anova; 3-summary & normality_test & multiple_comparisons]]

fitting.yb <- function(df, Col_x, Col_y, Col_group) {
  out <- data.frame()
  group <- as.character(unique(df[,Col_group]))
  for (i in c(1:length(Col_y))) {
    df1 <- data.frame(x = df[,Col_x],
                      y = df[,Col_y[i]],
                      group = df[,Col_group])
    df1$x <- as.numeric(as.character(df1$x))*100
    for (j in c(1:length(group))) {
      df2 <- df1[df1$group == group[j],]
      aa <- lm(y ~ x, data = df2) %>%
        summary()
      bb <- c(Col_y[i],group[j],
              round(aa$r.squared,3), 
              round(pf(aa$fstatistic[1],aa$fstatistic[2],aa$fstatistic[3],lower.tail=F),3))
      out <- rbind(out, bb)
    }
  }
  names(out) <- c("DV", "Treat", "R2", "P")
  out$sig[out$P < 0.05] <- "*"
  out$sig[out$P >= 0.05] <- ""
  return(out)
}
#直线回归：输入：df-数据表(数据框)、Col_x-自变量/x轴(1个列名)、Col_y-因变量/y轴(列名向量)、group-颜色(1个列名，4个颜色)
#输出：数据框-R^2、P、sig-P < 0.05

fitting.fig.yb <- function(df, Col_x, Col_y, Name_y, Col_group) {
  library(ggplot2)
  library(ggrepel)
  library(plyr)
  color <- c("#8dd3c7", "#fdb462", "#bebada", "#CC7D6D", "#e41a1c", "#377eb8")
  group <- unique(df[,Col_group])
  fig <- list()
  for (i in c(1:length(Col_y))) {
    df1 <- data.frame(x = df[,Col_x],
                      y = df[,Col_y[i]],
                      group = df[,Col_group])
    df1$x <- as.numeric(as.character(df1$x))*100
    fig[[i]] <- ggplot()+
      scale_color_manual(values = color)+
      labs(x = "", y = Name_y[i])+
      theme_bw()
    myp1 <- data.frame()
    for (j in c(1:length(group))) {
      df2 <- df1[df1$group == group[j],]
      aa <- lm(y ~ x, data = df2) %>%
        summary()
      bb <- c(aa$r.squared, pf(aa$fstatistic[1],aa$fstatistic[2],aa$fstatistic[3],lower.tail=F))
      if (bb[2] < 0.05) {
        fig[[i]] <- fig[[i]]+
          geom_smooth(data = df2, aes(x=x, y=y), method = lm, formula = y ~ x, se = F, 
                      color = color[j], linetype = 1, size = 1.5)
      } else {
        ret <- loess(y ~ x, data = df2, span = 0.75)
        newX=seq(25,100,0.1)
        ret_p <- data.frame(x=newX, y=predict(ret, newdata=data.frame(x=newX)))
        myborders <- newX[which(diff(diff(ret_p$y)>0)!=0L)+1L]
        fig[[i]] <- fig[[i]]+
          geom_line(data = ret_p, aes(x=x, y=y), color = color[j], linetype = 2, size = 1.5)
        if (length(myborders) > 0) {
          k=diff(diff(ret_p$y)>0L) 
          k=k[k!=0L]
          mystat=k[1]>0
          mystat0=mystat
          len=length(myborders)
          myborders <- (myborders[1:len-1]+myborders[2:len])/2
          myp <- sapply(as.data.frame(rbind(c(newX[1],myborders),
                                            c(myborders,newX[length(newX)]))),
                        function(k){mystat <<- !mystat;
                        optimize(f=function(x) {
                          predict(ret, newdata=data.frame(x))
                        },
                        maximum=mystat, interval=k)
                        }
          )
          myp=as.data.frame(t(myp))
          names(myp) <- c("x", "y")
          myp$x=as.numeric(myp$x)
          myp$y=as.numeric(myp$y)
          myp$type= 0;
          if(mystat0) {
            myp$type = myp$type + rep(c(1,0), length.out = nrow(myp))
          } else {
            myp$type = myp$type + rep(c(0,1), length.out = nrow(myp))
          }
          myp$xlabel <- paste(sprintf("%0.1f", myp$x),"%")
          myp1 <- rbind(myp1, myp)
        }
      }
    }
    fig[[i]] <- fig[[i]]+
      geom_point(data = myp1, aes(x, y), color = color[myp1$type + 5], size = 4, shape = 18)+
      geom_text_repel(data = myp1, aes(x, y, label = xlabel), 
                      color = color[myp1$type + 5], 
                      segment.color = color[myp1$type + 5], 
                      max.overlaps = 10000, fontface = "bold")+
      theme(legend.position = "top",
            axis.title.y = element_text(size = 14),
            axis.text.y = element_text(size = 12),
            axis.text.x = element_text(size = 10, face = "bold"),
            axis.title.x = element_blank())
  }
  return(fig)
}
#拟合曲线；输入：df-数据表(数据框)、Col_x-自变量/x轴(1个列名)、Col_y-因变量/y轴(列名向量)、Name_y-y轴名、group-颜色(1个列名，4个颜色)
#输出图：散点+拟合线（若直回归线p<0.05为直线实线回归线，若直回归线p>=0.05为loess回归虚线)

fig.yb <- function(list, y) {
  library(ggplot2)
  color <- c("#8dd3c7", "#fdb462", "#bebada", "#CC7D6D", "#e41a1c", "#377eb8")
  fig <- list()
  for (i in c(1:length(list[[1]]))) {
    fig_data <- list[[1]][[i]][[3]]
    fig_data$Treat <- paste(fig_data$T, fig_data$U, sep = "-")
    fig_data$Treat[fig_data$Treat == "0-0"] <- "CK"
    fig_data$Treat[fig_data$Treat == "1-0"] <- "T"
    fig_data$Treat[fig_data$Treat == "0-1"] <- "U"
    fig_data$Treat[fig_data$Treat == "1-1"] <- "T × U"
    fig_data$Treat <- factor(fig_data$Treat, levels = c("CK", "T", "U", "T × U"))
    fig[[i]] <- ggplot(fig_data)+
      geom_col(aes(x = Treat, y = mean, color = ID, fill = ID), width = 0.75, position = position_dodge())+
      geom_errorbar(aes(x = Treat, ymin = down, ymax = up, color = ID), 
                    width = 0.2, position = position_dodge(.75), size = 1)+
      geom_text(aes(x = Treat, y = up, label = HSD),
                position = position_dodge2(.75), vjust = -0.5, size = 3.5)+
      scale_color_manual(values = color)+
      scale_fill_manual(values = color)+
      scale_y_continuous(expand = c(0,0), limits = c(0,max(fig_data$up)*1.15))+
      labs(y = y[i])+
      theme_bw()+
      theme(legend.position = "none",
            axis.text.x = element_text(size = 12, face = "bold"),
            axis.title.x = element_blank(),
            axis.title.y = element_text(size = 14),
            axis.text.y = element_text(size = 12))
  }
  return(fig)
}

transpose.yb <- function(list, DV) {
  New <- list()
  for (i in c(1:length(list))) {
    for (j in c(1:length(list[[i]]))) {
      for (k in c(1:length(list[[i]][[j]]))) {
        dataf <- cbind(data.frame(DV = rep(DV[j], length(list[[i]][[j]][[k]][,1]))),
                       list[[i]][[j]][[k]])
        if (i == 1 & j == 1) {
          New[[k]] <- dataf
        } else {
          New[[k]] <- rbind(New[[k]], dataf)
        }
      }
    }
  }
  return(New)
}     #将多重列表转化为一层列表，用于输出excel

##################
###  数据分析  ###
##################
#二级数据计算
grow$TB <- grow$AB + grow$UB                                              #总生物量
grow$RS <- grow$UB / grow$AB                                              #根冠比

#soil$CN <- soil$SC / soil$SN                                              #碳氮比
#soil$CP <- soil$SC / soil$SP                                              #碳磷比
#soil$NP <- soil$SN / soil$SP                                              #氮磷比

result_O_grow <- ANOVA.yb(grow, c("H", "RL", "BS", "TB", "RS"), c("T", "U", "ID"))
result_O_leaf <- ANOVA.yb(leaf, c("leaf_a"), c("T", "U", "ID"))
result_O_ps <- ANOVA.yb(ps, c("A", "Ci", "Fv.Fm"), c("T", "U", "ID"))
#result_O_soil <- ANOVA.yb(soil, c("SC", "SN", "SP", "CN", "CP", "NP"), c("T", "U", "ID"))

grow_t <- Normal_t.yb(grow, c("H", "RL", "BS", "TB", "RS"), c("T", "U", "ID"))
leaf_t <- Normal_t.yb(leaf, c("leaf_a"), c("T", "U", "ID"))
ps_t <- Normal_t.yb(ps, c("A", "Ci", "Fv.Fm"), c("T", "U", "ID"))
#soil_t <- Normal_t.yb(soil, c("SC", "SN", "SP", "CN", "CP", "NP"), c("T", "U", "ID"))

result_N_grow <- ANOVA.yb(grow_t[[1]], c("H", "RL", "BS", "TB", "RS"), c("T", "U", "ID"))
result_N_leaf <- ANOVA.yb(leaf_t[[1]], c("leaf_a"), c("T", "U", "ID"))
result_N_ps <- ANOVA.yb(ps_t[[1]], c("A", "Ci", "Fv.Fm"), c("T", "U", "ID"))
#result_N_soil <- ANOVA.yb(soil_t[[1]], c("SC", "SN", "SP", "CN", "CP", "NP"), c("T", "U", "ID"))

result_grow <- list()
result_grow[[1]] <- list()
for (i in c(1:length(result_O_grow[[1]]))) {
  result_grow[[1]][[i]] <- list()
  result_grow[[1]][[i]][[1]] <- result_N_grow[[1]][[i]][[1]]
  result_grow[[1]][[i]][[2]] <- result_N_grow[[1]][[i]][[2]]
  result_grow[[1]][[i]][[3]] <- cbind(result_O_grow[[1]][[i]][[3]][,c(1:13)],
                                      result_N_grow[[1]][[i]][[3]][,c(14:21)])
}

result_leaf <- list()
result_leaf[[1]] <- list()
for (i in c(1:length(result_O_leaf[[1]]))) {
  result_leaf[[1]][[i]] <- list()
  result_leaf[[1]][[i]][[1]] <- result_N_leaf[[1]][[i]][[1]]
  result_leaf[[1]][[i]][[2]] <- result_N_leaf[[1]][[i]][[2]]
  result_leaf[[1]][[i]][[3]] <- cbind(result_O_leaf[[1]][[i]][[3]][,c(1:13)],
                                      result_N_leaf[[1]][[i]][[3]][,c(14:21)])
}

result_ps <- list()
result_ps[[1]] <- list()
for (i in c(1:length(result_O_ps[[1]]))) {
  result_ps[[1]][[i]] <- list()
  result_ps[[1]][[i]][[1]] <- result_N_ps[[1]][[i]][[1]]
  result_ps[[1]][[i]][[2]] <- result_N_ps[[1]][[i]][[2]]
  result_ps[[1]][[i]][[3]] <- cbind(result_O_ps[[1]][[i]][[3]][,c(1:13)],
                                    result_N_ps[[1]][[i]][[3]][,c(14:21)])
}

Table <- c(transpose.yb(result_grow, c("H", "RL", "BS", "TB", "RS")),
           transpose.yb(result_leaf, c("leaf_a")),
           transpose.yb(result_ps, c("A", "Ci", "Fv.Fm")))
Table[[1]] <- rbind(Table[[1]], Table[[4]], Table[[7]])
Table[[2]] <- rbind(Table[[2]], Table[[5]], Table[[8]])
Table[[3]] <- rbind(Table[[3]], Table[[6]], Table[[9]])
Table[[4]] <- rbind(grow_t[[2]], leaf_t[[2]], ps_t[[2]])
sheetname <- c("SD", "ANOVA", "MC", "Transformations")

if (file.exists(paste("Results ", Sys.Date(), ".xlsx", sep = "-"))) {
  file.remove(paste("Results ", Sys.Date(), ".xlsx", sep = "-"))
}
for (i in 1:length(sheetname)) {
  if (i == 1) {write.xlsx(Table[[i]], file = paste("Results ", Sys.Date(), ".xlsx", sep = "-"),
                          sheetName = sheetname[i])}
  else {write.xlsx(Table[[i]], file = paste("Results ", Sys.Date(), ".xlsx", sep = "-"),
                   sheetName = sheetname[i], append = TRUE)}
}

fit_grow <- fitting.yb(grow, "ID", c("H", "RL", "BS", "TB", "RS"), "Treat")
fit_leaf <- fitting.yb(leaf, "ID", c("leaf_a"), "Treat")
fit_ps <- fitting.yb(ps, "ID", c("A", "Ci", "Fv.Fm"), "Treat")

fig <- list()
color <- c("#8dd3c7", "#fdb462", "#bebada", "#CC7D6D", "#e41a1c", "#377eb8")
fig[[1]] <- c(fig.yb(result_grow, c("Height (cm)", "Root length (cm)", 
                                    "Stem diameter (mm)", "Total biomass (g)", "Root:shoot ratio")),
              fig.yb(result_leaf, c(expression("Leaf area (" * mm ^ 2 * ")"))),
              fig.yb(result_ps, c(expression(atop('Net photosynthetic rate', paste('(', μmol ~ m ^ -2 ~ s ^ -1, ')'))), 
                                  expression("Intercellular " * CO[2] * " (ppm)"),
                                  "Fv/Fm")))
Fig1_legend_data <- data.frame(Treat = c(1:4), mean = c(1:4), 
                               ID = as.factor(c("25%", "50%", "75%", "100%")))
Fig1_legend <- cowplot::get_legend(
  ggplot(Fig1_legend_data)+
    geom_col(aes(x = Treat, y = mean, color = ID, fill = ID), width = 0.75, position = position_dodge())+
    scale_color_manual(values = color, 
                       name = "Degree of invasion",
                       breaks = c("25%", "50%", "75%", "100%"),   
                       labels = c("25%", "50%", "75%", "100%"))+
    scale_fill_manual(values = color, 
                      name = "Degree of invasion",
                      breaks = c("25%", "50%", "75%", "100%"),   
                      labels = c("25%", "50%", "75%", "100%"))+
    theme(legend.position = "top",
          legend.title = element_text(size = 12, face = "bold"),
          legend.text = element_text(size = 12, face = "bold"))
)
Fig1 <- annotate_figure(
  ggarrange(fig[[1]][[1]], fig[[1]][[2]], fig[[1]][[3]],
            fig[[1]][[4]], fig[[1]][[5]], fig[[1]][[6]],
            fig[[1]][[7]], fig[[1]][[8]], fig[[1]][[9]], 
            align = "hv", labels = LETTERS, font.label = list(size = 16)),
  top = Fig1_legend)


fig[[2]] <- c(fitting.fig.yb(grow, "ID", 
                             c("H", "RL", "BS","TB", "RS"), 
                             c("Height (cm)", "Root length (cm)", 
                               "Stem diameter (mm)", "Total biomass (g)", "Root:shoot ratio"), 
                             "Treat"),
              fitting.fig.yb(leaf, "ID", 
                             c("leaf_a"), 
                             c(expression("Leaf area (" * mm ^ 2 * ")")), 
                             "Treat"),
              fitting.fig.yb(ps, "ID", 
                             c("A", "Ci", "Fv.Fm"), 
                             c(expression(atop('Net photosynthetic rate', paste('(', μmol ~ m ^ -2 ~ s ^ -1, ')'))), 
                               expression("Intercellular " * CO[2] *" (ppm)"),
                               "Fv/Fm"), 
                             "Treat"))

Fig2_x <- tableGrob("Degree of invasion", theme = ttheme_minimal(core=list(base_size = 16, fg_params=list(fontface="bold"))), cols = NULL)
Fig2_legend_data <- data.frame(x = rep(c(1:6),2), y = rep(c(1:6),2), 
                               type = as.factor(rep(c("CK", "T", "U", "T × U", 
                                    "Inflection point (+)", "Inflection point (-)"),2)))
Fig2_legend <- cowplot::get_legend(
  ggplot()+
    geom_line(data = Fig2_legend_data[c(1:4, 7:10),], 
              aes(x = x, y = y, color = type), size = 1.5)+
    #geom_point(data = Fig2_legend_data[5:6,], aes(x = x, y = y, color = type))+
    #geom_text(data = Fig2_legend_data[5:6,], aes(x = x, y = y, color = type, label = type))+
    scale_color_manual(values = color, 
                       breaks = c("CK", "T", "U", "T × U"),   
                       labels = c("CK", "T", "U", "T × U"))+
    theme_bw()+
    theme(legend.position = "top",
          legend.title=element_blank(),
          legend.text = element_text(size = 12, face = "bold"))
)

Fig2 <- annotate_figure(
  ggarrange(fig[[2]][[1]], fig[[2]][[2]], fig[[2]][[3]],
            fig[[2]][[4]], fig[[2]][[5]], fig[[2]][[6]],
            fig[[2]][[7]], fig[[2]][[8]], fig[[2]][[9]], 
            align = "hv", labels = LETTERS, font.label = list(size = 16)),
  bottom = Fig2_x, top = Fig2_legend)
  
Graphical_abstracts <- ggarrange(
  annotate_figure(fig[[1]][[4]]),
  annotate_figure(fig[[2]][[4]]),
  align = "hv", ncol = 2
)


###存图
#ggsave(file = "Fig1.svg", plot = Fig1, width = 16, height = 10,  bg = "white")
#ggsave(file = "Fig2.svg", plot = Fig2, width = 16, height = 10,  bg = "white")


#ggsave(file = "Fig1.tiff", plot = Fig1, width = 16, height = 10,  bg = "white")
#ggsave(file = "Fig2.tiff", plot = Fig2, width = 16, height = 10,  bg = "white")


ggsave(file = "Fig2-0315.png", plot = Fig1, width = 16, height = 10,  bg = "white")
ggsave(file = "Fig3-0315.png", plot = Fig2, width = 16, height = 10,  bg = "white")
ggsave(file = "Graphical_abstracts-0315.png", plot = Graphical_abstracts, width = 8, height = 3.7,  bg = "white")
